Event History Analysis III – Cox proportional hazard regression

Dagens program

• Cox-regression.
  • Opsummering på sidste gang uden tidsvarierende variable.
• Introduktion til tidsvarierende variable.
• Introduktion og diskussion af data i forlængelse af sidste lektion.
• Cox-regression med tidsvarierende variable.
• Øvelser.

Forventet udbytte

Opsummering

Spg.: Hvad er forskellen(e) på KM og Cox PH

Cox PH: fixed covariater

Hazarden for den \(i\)’ende person til tiden \(t\), er et produkt af to led i:

\[ h_{i}(t)=h_{0}(t) \times e^{\beta_{1} x_{1} + \dots + \beta_{k} x_{k}} \]

Cox PH: fixed covariater

Hazarden for den \(i\)’ende person til tiden \(t\), er et produkt af to led i:

\[ \require{color} h_{i}(t) = \;\;\; \textcolor{red}{h_{0}(t)} \;\;\; \times e^{\beta_{1} x_{1} + \dots + \beta_{k} x_{k}} \]

\(h_{0}(t)\) er baseline hazard ved tiden \(t\) for en referenceperson (en person med \(x=0\)). Baseline hazard specificeres ikke på forhånd, hvilket gør modellen semi-parametrisk. Bemærk, at der ikke er et separat intercept, fordi denne information er “absorberet” i \(h0(t)\). (Se også slides fra forrige lektion).

Cox PH: fixed covariater

Hazarden for den \(i\)’ende person til tiden \(t\), er et produkt af to led i:

\[ \require{color} h_{i}(t) = h_{0}(t) \times \;\;\; \textcolor{yellow}{e^{\beta_{1} x_{1} + \dots + \beta_{k} x_{k}}} \;\;\; \]

\(h_{0}(t)\) er baseline hazard ved tiden \(t\) for en referenceperson (en person med \(x=0\)). Baseline hazard specificeres ikke på forhånd, hvilket gør modellen semi-parametrisk. Bemærk, at der ikke er et separat intercept, fordi denne information er “absorberet” i \(h_{0}(t)\). (Se også slides fra forrige lektion).

\(e^{\beta_{k} x_{k}}\) angiver effekten af de \(k\) covariater. Effekten er konstant over tid og multiplicerer (på) baseline hazard for at give den samlede risiko for den enkelte person med specifikke værdier for covariaterne.

Cox PH: forudsætninger

Forudsætning for Cox ’s regressionsmodel er proportionalitet mellem hazard-raterne (den proportionale hazard model).

For to grupper gælder det at:

\[ h_{i}(t) = \theta h_{j}(t) \]

Derfor: ved sammenligning af hazard-raterne for to personer (\(i\) og \(j\)) på sammentidspunkt gælder det for hazard-ration at:

\[ \hat{HR} = \frac{h_{i}(t)}{h_{j}(t)} = \theta \] hvor \(\theta\) er en konstant.

Cox PH: Hazard-ration

Hazard-ratioen sammenligner risikoen for en begivenhed for en person uden for referencegruppen (\(i\)) med en person i referencegruppen (\(j\)), udtrykt som:

\[ \require{cancel} HR = \frac{h_{i}(t)}{h_{j}(t)} = \frac{ h_{0}(t) \times \text{exp}(\beta)}{h_{0}(t)} = \frac{ \cancel{h_{0}(t)} \times \text{exp}(\beta)}{\cancel{h_{0}(t)}} = \text{exp}(\beta) \]

Bemærk at \(h_{0}(t)\) ophæves, så HR kun afhænger af \(\text{exp}(\beta)\). Dette svarer konceptuelt til, hvordan Odds Ratio (OR) fungerer i logistisk regression.

• Hvis HR er \(=1\) har personen uden for referencegruppen samme risiko som personen i referencegruppen.

• Hvis HR er \(<1\) er risikoen lavere for personen uden for referencegruppen.

• Hvis HR er \(>1\) er isikoen højere for personen uden for referencegruppen.

Da HR er et forhold (en ratio), udtrykker den odds—ikke direkte sandsynligheder.

Cox PH: Hazard-ration og proportionalitet

library(tidyverse)

# Tid:
t <- seq(0, 10, by = 0.1)

# beta-værdien: effekten af covariaten, k:
beta <- 0.5

# dataframe:
df <- tibble(
  tid = t,
  h0 = 0.1 + 0.05 * sin(t), # Baseline hazard
  h1 = h0 * exp(beta),      # Hazard for gruppe i
  HR = h1 / h0              # Hazard-ratio mellem i og j
)

# Plot:
p1 <- ggplot(df, aes(x = tid)) +
  geom_line(aes(y = h0, color = "j"), size = 1.2) +
  geom_line(aes(y = h1, color = "i"), size = 1.2) +
  labs(x = "Time", y = "Hazard", title = "Proportionale Hazardfunktioner", color = "Gruppe") +
  theme_minimal()

p1

library(tidyverse)

# Tid:
t <- seq(0, 10, by = 0.1)

# Baseline hazard-funktion for referencegruppen, j:
h0 <- 0.1 + 0.05 * sin(t)

# Beta (effekten) + "gamma" (en tidsafhængig effekt der varierer over tid):
beta <- 0.5
gamma <- 0.2

# Hazard-funktion for gruppe i med tidsafhængig ikke-proportional effekt:
h2 <- h0 * exp(beta + gamma * t)

# Hazard-ratioen:
HR <- h2 / h0  # Dette svarer til exp(beta + gamma*t) ... altså, jeg har lavet en teoretisk funktion, der er ikke-proportionel over tid:

# Dataframe:
df <- tibble(
  tid = t,
  h0 = h0,
  h2 = h2,
  HR = HR
)

# Plot:
p1 <- ggplot(df, aes(x = tid)) +
  geom_line(aes(y = h0, color = "j"), size = 1.2) +
  geom_line(aes(y = h2, color = "i"), size = 1.2) +
  labs(x = "Tid", y = "Hazard", title = "Ikke-Proportionale Hazardfunktioner", color = "Gruppe") +
  theme_minimal()

p1

Cox PH: Hazard-ration og proportionalitet

library(tidyverse)

# Tid:
t <- seq(0, 10, by = 0.1)

# beta-værdien: effekten af covariaten, k:
beta <- 0.5

# dataframe:
df <- tibble(
  tid = t,
  h0 = 0.1 + 0.05 * sin(t), # Baseline hazard
  h1 = h0 * exp(beta),      # Hazard for gruppe i
  HR = h1 / h0              # Hazard-ratio mellem i og j
)

# Plot:
p1 <- ggplot(df, aes(x = tid)) +
  geom_line(aes(y = h0, color = "j"), size = 1.2) +
  geom_line(aes(y = h1, color = "i"), size = 1.2) +
  labs(x = "Tid", y = "Hazard", title = "Proportionale Hazardfunktioner", color = "Gruppe") +
  theme_minimal()

p1

# Plot:
p2 <- ggplot(df, aes(x = tid, y = HR)) +
  geom_line(color = "green", size = 1.2) +
  labs(x = "Time", y = "Hazard Ratio", title = "Hazard Ratio (Konstant over tid)") +
  theme_minimal()

p2

library(tidyverse)

# Tid:
t <- seq(0, 10, by = 0.1)

# Beta (effekten) + "gamma" (en tidsafhængig effekt der varierer over tid):
beta <- 0.5
gamma <- 0.2

# Dataframe:
df <- tibble(
  tid = t,
  h0 = 0.1 + 0.05 * sin(t),
  h2 = h0 * exp(beta + gamma * t), # Hazard-funktion for gruppe i med tidsafhængig ikke-proportional effekt
  HR = h2 / h0  # Dette svarer til exp(beta + gamma*t) ... altså, jeg har lavet en teoretisk funktion, der er ikke-proportionel over tid
)

# Plot:
p1 <- ggplot(df, aes(x = tid)) +
  geom_line(aes(y = h0, color = "j"), size = 1.2) +
  geom_line(aes(y = h2, color = "i"), size = 1.2) +
  labs(x = "Tid", y = "Hazard", title = "Ikke-Proportionale Hazardfunktioner", color = "Gruppe") +
  theme_minimal()

p1

# Plot:
p2 <- ggplot(df, aes(x = tid, y = HR)) +
  geom_line(color = "green", size = 1.2) +
  labs(x = "Tid", y = "Hazard Ratio", title = "Hazard Ratio (Ikke-konstant over tid)") +
  theme_minimal()

p2

Tied data

Tied data opstår ofte i overlevelsesanalyser, fordi vi typisk måler tid i diskrete enheder. Det betyder, at flere hændelser kan ske i samme tidsinterval. Denne situation bryder med antagelsen om, at alle observationer er uafhængige, og kan derfor føre til bias i estimater og standardfejl, som påvirker validiteten af vores statistiske konklusioner.

Vi kan løse dette problem ved hjælp af tre forskellige standard-metoder:

• Breslow: Velegnet til “små” datasæt med få ties.

• Efron: Klarer bedre mange ties.

• Exact: Den mest nøjagtige metode, men også den mest beregningstunge (computerkrævende).

Tied data: Efron

Når flere personer oplever hændelsen i samme tidsinterval, kaldes det for tied data. Efron-metoden håndterer dette ved at tage højde for alle mulige rækkefølger (rank orderings) af de hændelser, der sker på samme tid. På den måde justeres, hvor meget de “tied” hændelser bidrager til den samlede likelihood-funktion når vi estimerer vores model.

Vi udtrykke med denne metode udtrykke hazardfunktion ved en distinkt tid \(t_{(i)}\)​ med tied events som:

\[ h_{a}(t_{(i)}) = \frac{\sum_{j=1}^{d_i} \delta_{ij}}{\sum_{j \in R(t_{(i)})} \sum_{k=1}^{d_k} \exp(X_j'\beta) \times \delta_{ik}} \]

hvor

• \(h_{a}(t_{(i)})\) er den kontrollerede/justerede hazard funktion ved \(t_{(i)}\),
• \(t_{(i)}\) er den distinkte tid, hvor hændelser opstå,
• \(j\) er indekset for individer som oplevede begivenheden ved \(t_{(i)}\),
• \(d_{i}\) er antallet af individer der oplevede begivenheden ved tid \(t_{(i)}\),
• \(R(t_{(i)})\) er risikosættet før \(t_{(i)}\),
• \(k\) er indekset for individer \(i\) risikosættet, \(R(t_{(i)})\), før \(t_{(i)}\),
• \(\delta_{ij}\) er en dummy, der indikerer hvis individ \(j\) oplevede begivenheden ved \(t_{(i)}\) (1) eller ej (0),
• \(X_{j}\) er en vector af covariater for individ \(j\),
• \(\beta\) er en vector af coeficcienter fra Cox PH modellen.

Tied data: Efron

Metoden tager højde for alle mulige rækkefølger af disse forbundne begivenheder inden for et givent tidsinterval. For eksempel:

• 3 personer, \(A, B, C\), oplevede begivenheden ved \(t=x\).
• Da vi arbejder med diskrete tidsenheder er der en række af mulige “korrekte” rækkefølger som \(A\), \(B\) og \(C\) oplevede begivenheden:
  • \(A\) først, så \(B\), så \(C\)
  • \(B\), så \(C\), så \(A\),
  • osv. …

På den måde “glatter”/smoother vi estimatet af hazardfunktionen og kontrollerer for den usikkerhed, der opstår, når begivenheder registreres som værende samtidige og observationerne potentielt er afhængige.

Tied data: Efron

Det vil sige, at vi overvejer alle mulige måder, de individer, der oplever begivenheden samtidigt, kan rangeres, hvilket justerer for den usikkerhed, der opstår ved tied data. Derfor får vi en mere præcis likelihood funktion, udtrykt som:

\[ L(\beta) = \prod_{i=1}^m \frac{\exp\left(\sum_{j \in R(t_{(i)})} X_j'\beta \right)^{d_i}}{\left[\sum_{j \in R(t_{(i)})} \exp(X_j'\beta) \right]^{d_i}} \]

hvor produktet over \(i\) går fra 1 to \(m\). \(m\) er det samlede antal diskrete begivenhedstindspunkter.

Fordi vi kun tager de individer, \(j\), der oplevede begivenheden i betragtning og ikke specificerer baseline hazard funktionen er det formelt en partial likelihood funktion.

Tied data: Exact

På mange måder ligner Exact metode Efron metoder. Den grundlæggende forskel er præcisionen i estimater (og den direkte relaterede computerkraft det kræver at være mere præcis). Med andre ord er \(h(t)\) ikke bare “kontrollere”/adjusted som med Efron metoden, men præcis.

Hvor Efron kontrollerede for mulige permutationer af forbundne begivenheder indenfor et sæt af forbundne begivenheder, tager “exact” metode alle mulige permutationer i dataen (og ikke bare sættet) i betragtning og er derfor ikke afhængig af et “smoothed” estimat men giver et præcist estimat.

Cox PH med tidsvariende variable

For at udvide Cox modellen med tidsvarierende covariate er der kun ganske få ændringer i formlen I er ved at være bekendt med:

\[ h_{i}(t)=h_{0}(t) \times e^{\beta_{1} x_{1} + \beta_{2} x_{2}\textcolor{yellow}{(t)} + \dots + \beta_{k} x_{k}} \]

Den eneste forskel er at vi tilføjer \(t\) for de tidsvarierende covariater.

Her er \(x_{2}(t)\) en tidsvarierende variabel, mens de øvrige \(x\)-variable er faste, og skal læses som at hazarden ved tid \(t\) er bestemt af tidsinvariante karakteristika, \(x_{1}\), og \(x_{2}\) ved tid \(t\).

Cox PH med tidsvariende variable

Nu gælder det at:

\[ \hat{HR} = \frac{h_{i}(t)}{h_{j}(t)} = \frac{h_{0}(t) \times \exp(\beta x_{i}\textcolor{yellow}{(t)})}{h_{0}(t) \times \exp(\beta x_{j}\textcolor{yellow}{(t)})} \]

Her ophæves \(h_{0}(t)\) fortsat, men fordi \(x_{i}(t)\) og \(x_{j}(t)\) kan ændre sig over tid, bliver HR nu afhængig af tidspunktet:

Med andre ord, hazardratioen bliver afhængig af tidspunktet og kan variere over tid, hvis de tidsvarierende variable for de to personer ændrer sig:

\[ \hat{HR} = \exp\left(\beta \left[x_i(t) - x_j(t)\right]\right) \]

Det betyder, at hazard-ratioen i denne situation bliver tidsafhængig, fordi den afhænger af, hvordan den tidsvarierende variabel adskiller sig mellem personerne på et givet tidspunkt.

Cox PH med tidsvariende variable: proportionalitet?

library(tidyverse)

# Tid
t <- seq(0, 10, by = 0.1)

# Beta: effekten af covariaten
beta <- 0.5

# Tidsvarierende covariater for de to grupper:
# Gruppe i: x_i(t) = 1 + 0.1*t
# Gruppe j: x_j(t) = 1 + 0.05*t
x_i <- 1 + 0.1 * t
x_j <- 1 + 0.05 * t

# Beregn hazard for hver gruppe:
h_i <- h0 * exp(beta * x_i)  # Hazard for gruppe i
h_j <- h0 * exp(beta * x_j)  # Hazard for gruppe j

# Opret en tibble med data
df <- tibble(
  tid = t,
  h0 = 0.1 + 0.05 * sin(t),
  `i` = h_i,
  `j` = h_j,
  HR = exp(beta * (x_i - x_j))
)

# Omform data til long format til ggplot
df_long <- df %>% pivot_longer(cols = c("i", "j"),
                               names_to = "Gruppe",
                               values_to = "Hazard")

# Plot:
p1 <- ggplot(df_long, aes(x = tid, y = Hazard, color = Gruppe)) +
  geom_line(size = 1.2) +
  labs(title = "Hazardfunktioner over tid for to grupper",
       x = "Tid",
       y = "Hazard") +
  theme_minimal()

p1

library(tidyverse)

# Tid
t <- seq(0, 10, by = 0.1)

# Beta: effekten af covariaten
beta <- 0.5

# Tidsvarierende covariater for de to grupper:
# Gruppe i: x_i(t) = 1 + 0.1*t
# Gruppe j: x_j(t) = 1 + 0.05*t
x_i <- 1 + 0.1 * t
x_j <- 1 + 0.05 * t

# Beregn hazard for hver gruppe:
h_i <- h0 * exp(beta * x_i)  # Hazard for gruppe i
h_j <- h0 * exp(beta * x_j)  # Hazard for gruppe j

# Opret en tibble med data
df <- tibble(
  tid = t,
  h0 = 0.1 + 0.05 * sin(t),
  `i` = h_i,
  `j` = h_j,
  HR = exp(beta * (x_i - x_j))
)

# Omform data til long format til ggplot
df_long <- df %>% pivot_longer(cols = c("i", "j"),
                               names_to = "Gruppe",
                               values_to = "Hazard")

# Plot: 
p2 <- ggplot(df, aes(x = tid, y = HR)) +
  geom_line(color = "green", size = 1.2) +
  labs(
    title = "Tidsvarierende Hazard Ratio",
    x = "Tid",
    y = "Hazard Ratio, HR(t)"
  ) +
  theme_minimal()

p2

Cox PH med tidsvariende variable: Data

flowchart TD

    X([<font size=6><font color=white> Afhængig <br/> variabel]) 

    A([<font size=6><font color=black> Forklarende <br/> variabel]) 
    -...-> B([<font size=6><font color=white>Tidsvariaende])
    
    A 
    ---> C([<font size=6><font color=black>Konstant])
    
    B
    --> D([<font size=6><font color=black> Variable der løber <br/> med tiden])
    
    B
    --> E([<font size=6><font color=black>En variabel med forskellige <br/> værdier til forskellige <br/> tidspunkter])
    
    style A fill:#f45c2c,stroke:#111,stroke-width:1px
    style B fill:#333,stroke:#111,stroke-width:1px
    style C fill:#d8897b,stroke:#111,stroke-width:1px
    style D fill:#9f6671,stroke:#111,stroke-width:1px
    style E fill:#5e4f6d,stroke:#111,stroke-width:1px
    

Cox PH med tidsvariende variable: Data

flowchart TD
    B([<font size=6><font color=white>Tidsvariaende])
    --> D([<font size=6><font color=black> Variable der løber <br/> med tiden])
    
    B
    --> E([<font size=6><font color=black>En variabel med forskellige <br/> værdier til forskellige <br/> tidspunkter])
    
    D 
    --- F([For eksempel alder])
    
    D 
    -...- G([<font size=6><font color=black>Kan altid anvendes bestemmes <br/> ved start- og slutoplysninger])
    
    E
    --- H([For eksempel arb. eller udd. status])
    
    E
    -...- I([<font size=6><font color=black>Skal bestemmes for <br/> hvert tidspunkt])
    
    style B fill:#d8897b,stroke:#111,stroke-width:1px
    style D fill:#9f6671,stroke:#111,stroke-width:1px
    style E fill:#5e4f6d,stroke:#111,stroke-width:1px
    style G fill:#b7bd50,stroke:#111,stroke-width:1px
    style I fill:#b7bd50,stroke:#111,stroke-width:1px

Cox PH med tidsvariende variable: Data

For at coxph() funktionen kan tage tidsvarierende informationer er det (desværre) ikke nok med periode_start og periode_slut. For hvert tidspunkt og person skal vi have start og slut information.

Cox PH med tidsvariende variable: Data

For at coxph() funktionen kan tage tidsvarierende informationer er det (desværre) ikke nok med periode_start og periode_slut. For hvert tidspunkt og person skal vi have start og slut information.

Cox PH med tidsvariende variable: Data

For at coxph() funktionen kan tage tidsvarierende informationer er det (desværre) ikke nok med periode_start og periode_slut. For hvert tidspunkt og person skal vi have start og slut information.

library(tidyverse)
library(readxl)

df.rå <- read_excel("eksempel_data.xlsx")

df.lagged <- 
  df.rå %>% 
  mutate(person_start = case_when(lag(id) == id ~ 0, TRUE ~ 1),
         person_slut = case_when(lead(id) == id ~ 0, TRUE ~ 1)) %>% 
  group_by(id) %>% 
  mutate(begivenhed = if_else(lag(børn) == 0 & børn == 1 & person_start == 0, 1, 0),
         periode_slut = case_when(begivenhed == 1 | børn == 0 & person_slut == 1 ~ 1, TRUE ~ 0)) %>% 
  slice(1:which.max(periode_slut == 1)) %>% 
  mutate(start = 1:n(),
         slut = start + 1) %>% 
  ungroup() 
  

Cox PH med tidsvariende variable: Data

Disse to variable gør det muligt at udvide coxph() funktionen med:

cox.model2 <- coxph(Surv(start, slut, event) ~ covariater, method = "efron", data = df)

Dette kalder vi også for et Counting process argument.

Her angives start og sluttiden for det enkelte interval, samt event (1=hændelse og 0=censurering). Intervallets størrelse afgøres af hvor hyppigt de tidsvarierende variable kan skifte værdi. Værdien skal være konstant inden for det enkelte interval. Denne lille ændring gør det muligt at håndtere tidsvarierende variable i cox-regressionen

Dagens øvelse: data

load("//Volumes//oecon-forloeb24//Råfiler//Forlob_14032023.rda")

Dagens øvelse: data

load("//Volumes//oecon-forloeb24//Råfiler//Forlob_14032023.rda")

Oprindeligt: - 1 % af alle individer fundet i IDAP igennem perioden 1980-2010

Reduktion: - Data indeholder personer, der er fyldt 10 år mellem 1980-1990 - Personer, der får børn før de er fyldt 18 år er taget ud af data (56 personer) - Kun personer, der er til stede i DK, det år de bliver 18 år. - Årlige observationer af disse personer